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ABSTRACT 

We investigate the characteristic radiative efficiency e, Eddington ratio A, and duty cycle Po of high-redshift 
active galactic nuclei (AGN), drawing on measurements of the AGN luminosity function at z = 3 — 6 and, 
especially, on recent measurements of quasar clustering at z = 3 — 4.5 from the Sloan Digital Sky Survey. The 
free parameters of our models are e. A, and the normalization, scatter, and redshift evolution of the relation 
between black hole mass Mbh and halo virial velocity Vvi,. We compute the luminosity function from the 
implied growth of the black hole mass function and the quasar correlation length from the bias of the host 
halos. We test our adopted formulae for the halo mass function and halo bias against measurements from 
the large N-body simulation developed by the MICE collaboration. The strong clustering of AGNs observed 
at z = 3 and, especially, at z = 4 implies that massive black holes reside in rare, massive dark matter halos. 
Reproducing the observed luminosity function then requires high efficiency e and/or low Eddington ratio A, 
with a lower limit (based on 2a agreement with the measured z = 4 correlation length) e > 0.7A/(1 + 0.7A), 
implying e > 0.17 for A > 0.25. Successful models predict high duty cycles, Pq ~ 0.2,0.5, and 0.9 at z = 3.1,4.5 
and 6, respectively, and they require that the fraction of halo baryons locked in the central black hole is much 
larger than the locally observed value. The rapid drop in the abundance of the massive and rare host halos at 
z > 7 implies a proportionally rapid decline in the number density of luminous quasars, much stronger than 
simple extrapolations of the z = 3 — 6 luminosity function would predict. For example, our most successful 
model predicts that the highest redshift quasar in the sky with true bolometric luminosity L > 10"^^^ erg s^' 
should be at z ^ 7.5, and that all quasars with higher apparent luminosities would have to be magnified by 
lensing. 

Subject headings: cosmology: theory - black hole: formation - galaxies: evolution - quasars: general 



1. INTRODUCTION 

The masses of the central black holes in local galaxies 
are correlated with the luminosities, stellar and dynamical 
masses, and velocity dispersions of the galaxies in which they 
reside (e.g., Magorrian et al. 1998; Ferrarese & Merritt 2000; 
Gebhardt et al. 2000; McLure & Dunlop 2002; Marconi 
& Hunt 2003; Hai'ing & Rix 2004; Ferrarese & Ford 2005; 
Greene & Ho 2006; Graham 2007; Hopkins et al. 2007b; 
Graham 2008). The Mbh -c relation, together with the ob- 
served correlation between outer circular velocity and cen- 
tral velocity dispersion measured by several groups (e.g., Fer- 
rarese 2002; Baes et al. 2003; Pizzella et al. 2005; Buyle et 
al. 2006), implies a mean correlation between black hole mass 
and the mass or virial velocity of the host galaxy's dark matter 
halo, although with a possibly large scatter (e.g.. Ho 2007 a,b). 
Recent observational studies have attempted to constrain the 
evolution of the black holes and their host galaxies, by mea- 
suring the Mbh -c relation at < z < 3, finding only tentative 
evidence for larger black holes at fixed velocity dispersion 
or stellar mass (e.g., McLure et al. 2006; Peng et al. 2006; 
Shields et al. 2006; Treu et al. 2007; Shankar et al. 2009b) 
with respect to what is observed locally. However, such an 
evolution is difficult to detect given the limited sampling and 
bias effects involving these measurements (e.g., De Zotti et al. 
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2006; Lauer et al. 2007; Ho 2007a). Probing this evolution 
becomes even more difficult at z > 2 because luminous AGNs 
substantially outshine their hosts. 

Another way to probe the evolution of black holes and their 
host galaxies comes from clustering. Since more massive ha- 
los exhibit stronger clustering bias (Kaiser 1984; Mo & White 
1996), the clustering of quasars provides an indirect diagnos- 
tic of the masses of halos in which they reside (Haehnelt et 
al. 1998; Haiman & Hui 2001; Martini & Weinberg 2001; 
Wyithe & Loeb 2005), which in turn can provide informa- 
tion on black hole space densities, on duty cycles and life- 
times, and, indirectly, on the physical mechanisms of black 
hole feeding. Measuring the clustering as a function of red- 
shift and quasar luminosity probes the relation between AGN 
luminosity and host halo mass, thus constraining the distri- 
butions of Eddington ratios and radiative efficiencies which 
govern the accretion of black holes at different epochs and in 
different environments. The strong clustering of quasars at 
z > 3 recently measured by Shen et al. (2007; hereafter S07) 
in the Sloan Digital Sky Survey (SDSS; York et al. 2000) 
quasar catalog (Richards et al. 2002; Schneider et al. 2007) 
implies that the massive black holes powering these quasars 
reside in massive, highly biased halos. 

The classical modeling of quasar clustering by Haiman 
& Hui (2001) and Martini & Weinberg (2001) assumes a 
mean value for the duty cycle and derives the relation be- 
tween quasar luminosity and host halo mass by monotonically 
matching their cumulative distribution functions. White, Mar- 
tini & Cohn (2008; hereafter WMC) have appHed this method 
to the S07 measurements, concluding that the strong cluster- 
ing measured at z ^ 4 can be understood only if quasar duty 
cycles are high and the intrinsic scatter in the luminosity-halo 
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relation is small. In this paper, we take a further step by jointly 
considering the evolution of the black hole-halo relation and 
the black hole mass function, as constrained by the observed 
AGN luminosity function and clustering. We examine con- 
straints on the host halos, duty cycles, radiative efficiencies, 
and mean Eddington ratios of massive black holes at z > 3, 
imposed by the clustering measurements of S07 and by a va- 
riety of measurements of the quasar luminosity function at 
3 < z < 6 (e.g., Kennefick et al. 1994; Pei 1995; Fan et al. 
2001, 2004; Barger et al. 2003; Wolf et al. 2003; Hunt et al. 
2004; Barger & Cowie 2005; La Franca et al. 2005; Nandra et 
al. 2005; Cool et al. 2006; Richards et al. 2006a; Bongiorno 
et al. 2007; Fontanot et al. 2007; Shankar & Mathur 2007; 
Silverman et al. 2008; Shankar et al. 2010a,b). 

Our method of incorporating luminosity function con- 
straints is simple. We assume the existence of a relation be- 
tween black hole mass Mbh and halo virial velocity Vvh at 
high redshift and assume that the slope of this relation is the 
same as observed locally, but leave its normalization, redshift 
evolution and scatter as adjustable parameters. Since the halo 
mass function is predicted from theory at every redshift, the 
evolution of the black hole mass function follows once the 
Mbh -Kir relation is specified. This growth of black holes is 
then used to predict the AGN luminosity function, in terms of 
the assumed radiative efficiency e/(l — e) = L/Mbhc^ and Ed- 
dington ratio A = L/LEdd of black hole accretion, which can 
be compared to observations. This method inverts the "con- 
tinuity equation" approach to quasar modeling, in which one 
uses the observed luminosity function to compute the implied 
growth of the black hole mass function (e.g., Cavaliere et al. 
1971; Small & Blandford 1992; Yu & Tremaine 2002; Steed 
& Weinberg 2003; Mai'coni et al. 2004; Shankai- et al. 2004; 
Yu & Lu 2004; Shankar, Weinberg, & Miralda-Escude 2009a, 
hereafter SWM; Shankar 2009). 

We make no specific hypothesis about the mechanisms that 
trigger high-redshift quasar activity. Our model simply as- 
sumes that a relation between Mbh and Vvii- exists and that 
it is maintained by mass accretion that produces luminous 
quasar activity, assuming no significant time delay between 
the two. As detailed below, simultaneously matching the ob- 
served luminosity function and the S07 clustering measure- 
ments, especially their z = 4 correlation length, is in general 
quite difficult. Moderately successful models must share the 
common requisites of having low intrinsic scatter in the Mbh 
- Vvir relation and a high value of the ratio e/A. Although 
these findings are affected by the model adopted to compute 
the halo bias factor, we will show that they do not otherwise 
depend on the details of our modeling and can be understood 
in simple, general terms. 

The mass function and clustering bias of rare, massive ha- 
los at high redshift are crucial inputs for our modeling. We 
therefore test existing analytic formulae for these quantities 
against measurements from the large N-body simulation de- 
veloped by the MICE collaboration, which uses 10^ particles 
to model a comoving volume of 768/i^' Mpc on a side. 

Throughout this paper the following cosmological param- 
eters have been used, consistent with the best-fit model to 
WMAP5 data (Spergel et al. 2007): n„, = 0.25, fit = 0.75, 
(78=0.8, n = 0.95, /! = //o/100kms-iMpc-i =0.7, fib = 
0.044. 

2. MODEL 

2.1. AGN BIAS AND LUMINOSITY FUNCTION 



In the local universe, the masses of black holes are tightly 
correlated with the velocity dispersion a of their parent bulges 
(e.g., Fen-arese & Memtt 2000; Gebhardt et al. 2000; 
Tremaine et al. 2002; Shankar & Ferrarese 2009). This re- 
lation has been recently re-calibrated by Tundo et al. (2007) 
as 

logf^)=8.21+3.831og(^,;,;^V (1) 
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where we denote the average black hole mass at a fixed <t as 
Mbh- The bulge velocity dispersions are in turn correlated 
with large scale circular velocities (Ferrarese 2002; see also 
Baes et al. 2003 and Pizzella et al. 2005): 

log V, = (0.84 ± 0.09) logCT + (0.55 ±0.19), (2) 

with <T and Vc measured in kms^'. For a flat rotation curve, 
the disk circular velocity is equal to the halo virial velocity 
Vvii . Departures from isothermal halo profiles, gravity of the 
stellar component, and adiabatic contraction of the inner halo 
can alter the ratio Vvir /Vc, but the two quantities should re- 
main well correlated nonetheless (e. g.. Mo et al. 1998; Mo & 
Mao 2004 and references therein). Thus, the correlations ([U 
and (|2]l imply a correlation between black hole mass and halo 
virial velocity, although we should expect the Mbh -Kh rela- 
tion to have a larger scatter than the observed Mbh -o' relation 
(e.g.. Ho 2007b). 

As mentioned in §[1] the models we shall construct assume 
that black holes at z > 3 lie on an Mbh -Kir relation of similar 
form. We parameterize this relation as 
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which corresponds to equations ^ and dU with Vvir replacing 
Vc- We define a as the normalization of the Mbh -Kir rela- 
tion at z = 3.1, which corresponds to the mean redshift in the 
lower subsample of S07. The factor a allows both for an off- 
set between the z = 3 . 1 and z = relations and for a ratio V^r 
/Vc ^ I at z = 0. For example, typical disk galaxy models 
(e.g.. Mo et al. 1998; Seljak 2002; Dutton et al. 2007; Gnedin 
et al. 2007) have K/Kir « 1.4 - 1.8 at z = 0, which would 
imply normalizations a « 5 — 15 for equation (O at z = be- 
cause of the steep power of velocity. Note that none of our 
results depend on the z = normalization of the Mbh ~ cr re- 
lation because we use only high redshift data in this paper In 
addition, we allow redshift evolution in the Mbh -Kir relation 
at z > 3. 1 through the index 7. 

The relation between the halo virial velocity K-ir and the 
halo virial mass M is 
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where the mean density contrast (relative to critical) within 
the virial radius R^i, is A = ISn^ + 82c/ - 39d^, with 
d = n,„{z) - 1, and n,„{z) = a„(l +z)V [^n,{l +zy + ^Ja] 
(Bryan and Norman 1998; Barkana & Loeb 2001). In terms 
of halo mass, equation ^ corresponds to 
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We assume the presence of a scatter about this mean relation, 
with a log-normal distribution and a dispersion S in the loga- 
rithm of black hole mass at fixed Vvir . 

Given the theoretically known halo mass function, we com- 
pute the black hole mass function via the convolution 
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(logMBH[M,z] - logMBH)^ 
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with 
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where E is the log-normal scatter in Mbh at fixed halo mass, 
X =Mbh or M, and ns{x,z)dx is the comoving number den- 
sity of black holes/halos (for subscript s — BH or s = h) in 
the mass range x ^ x + dx, in units of Mpc^-' for Hq = 
70kms^' Mpc^'. The units of "I>.s are comoving Mpc^^ per 
decade of mass. We convert to these units in order to compare 
with the data on the AGN luminosity function. 

The quasar luminosity function $(L,z), expressed in the 
same units as is modeled according to a simple pre- 

scription where black holes can be in only two possible states: 
active or inactive. All black holes that are active accrete with 
a single value of the radiative efficiency, e, and of the Edding- 
ton ratio, A = L/LEdd, where L is the bolometric luminosity 
and 



LEdd = 1-26X 10^** ergs"' 
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is the Eddington luminosity (Eddington 1922). The growth 
rate of an active black hole of mass Mbh is Mbh =Mbh / fef, 
where the e-folding time is (Salpeter 1964) 



fef = 4x 10**( ^ ) yr, 



(9) 



where f — e/ {I — e), and the radiative efficiency is e = L(l — 
e)/[MBHC^]- (Radiative efficiency e is conventionally defined 
with respect to the mass inflow rate M, and the black hole 
mass growth rate Mbh is smaller by a factor 1 — e because of 
radiative losses). 

Once the parameters a and 7 of the Mbh -Vvir relation are 
specified, the growth of $bh(Mbh,z) is determined by the 
(theoretically calculable) evolution of the halo mass function 
n/,(M,z). We compute the AGN luminosity function assuming 
that this growth is produced by accretion with radiative effi- 
ciency e and Eddington ratio A. This method inverts a long- 
standing approach to modeling AGN and black hole evolution 
in which one calculates the growth of the black hole mass 
function implied by the observed luminosity function using a 
"continuity equation". 
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(see, e.g., Cavahere et al. 1971; Small & Blandford 1992; 
Yu & Tremaine 2002; Mai'coni et al. 2004; SWM). Here we 
ignore the impact of black hole mergers in the evolution of 
the black hole mass function, because the black hole mass 
growth via mergers is relatively small, as we show in detail in 
the Appendix. 



Knowing $bh(Mbh,z), we can invert equation (fTOl i to ob- 
tain the luminosity function 
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In practice, we integrate equation ( fTTT ) up to black hole masses 
of logMBH/M0 = 11. Equation ( fTTT i assumes a strictly mono- 
tonic, scatter-free relation between AGN luminosity and black 
hole mass. Therefore in our models the only source of scatter 
between AGN luminosity and halo mass is the scatter in the 
Mbh -Kir relation. However, provided that the L-Vvir scatter 
is fairly small (as we find it must be to explain the observed 
clustering), we expect that it makes little difference whether 
it arises from scatter in Mbh or scatter in A. 

The average growth rate of all black holes (active and inac- 
tive) of mass Mbh is (Mbh) = A)MBH/fef, where the duty cycle 
P() is the probability that a black hole is in the active state. In 
models with a single value of A, the duty cycle is simply the 
ratio of the luminosity function to the mass function. 
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«'bh(Mbh,z) 

where XI A physically consistent model must have 

P() < 1 for all Mbh and z, and can be directly computed from 
equations |6] and [TT| 

In addition to the AGN luminosity function, we test our 
models against the clustering measurements of S07, specifi- 
cally their reported values of the AGN correlation length tq. 
We calculate these correlation lengths from the condition 



Piz)D^iz)arQ) = l: 



(13) 



where ^(ro) is the Fourier transform of the linear power spec- 
trum, D{z) is the linear growth factor of perturbations, and 
b{z) is the mean clustering bias of AGN shining above a lu- 
minosity threshold L^in at redshift z, given by (Haiman & Hui 
2001) 
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The minimum luminosity Lmin (z) in equation ( fT4] i is a bolo- 
metric quantity, while the S07 bias is measured above a red- 
shift dependent, A'-corrected M, magnitude. To convert from 
magnitudes to bolometric luminosities, we first convert to B 
magnitudes assuming Mb — Mi{z — 2) +0.804 (Richards et 
al. 2006b), and then adopt an average bolometric correction 
of Cb — 10.4, with L = CbLbVb^ where vb is the frequency at 
the center of the B-band (at wavelength 4400 A). Because our 
models assume a single Eddington ratio A, b{L,z) is just equal 
to the bias Z7(Mbh,z) of black holes of mass Mbh = i/^ A. The 
latter is computed from the b{M,z) of halos of mass M using 
the model relation between Mbh and M(z) (equation 15] ). In- 
cluding the log-normal scatter of width E, the black hole bias 
is 

b{Msu,z) - [$BH(MBH,z)]"'y" b{M,z)MM,z) x 



(2^E2)-i/2exp 



(logMBH[M,z] -logMfin) 



2E2 



dlogHlS) 



We discuss our choice of b{M,z) in § 12.21 
In summary, the free parameters of our model are: 
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• the normalization constant a in the Mbh -Kir relation, 

• the parameter 7 which regulates the redshift evolution 
[ ( 1 + 2) /4 . 1 ] 7 of this relation, 

• the mean Eddington ratio A of active black holes, 

• the log-normal scatter S in Mbh at fixed Vvir , 

• the radiative efficiency e of black hole accretion. 

The predicted bias in these models is completely independent 
of the assumed radiative efficiency (when other parameters 
are held fixed), since the efficiency does not affect the relation 
between luminosity and halo mass. 
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Fig. 1. — Halo ma.ss function: compari.son between N-body measurements 
and analytic fits. Upper panel : Measured mass functions for FoF hales iden- 
tified at redshifts 3 (filled circles), 4 (filled sc/uares), 5 (open triangles) and 6 
(open circles) in our simulation are in better agreement with Sheth & Tormen 
(1999) (discontinuous lines) rather than the Jenkins et al. (2001) (solid lines) 
analytic fits. Fits are displayed for the same redshifts as measurements in the 
simulation: from z = 3 (upper lines) to z = 6 (bottom lines). Lower panel 
: Fractional deviations of the measured mass functions (symbols as in upper 
panel) and the Sheth & Tormen predictions (discontinuous lines as above) 
with respect to the Jenkins et al. fit (solid line). 



2.2. MASS FUNCTION AND HALO BIAS 

The high redshift quasars used by S07 and our study are be- 
lieved to reside in very rare halos with M ~ IO'^-i^/j-IMq 
at redshifts z = 3 to 6. While extensive work has been done 
to determine the abundances and clustering of halos at z < 3, 
testing the accuracy of simple analytic formulae against pre- 
dictions from cosmological numerical simulations of struc- 
ture formation, this work has not been extended to the high- 
redshift (z ^ 3 — 6), rare massive halos we are interested in 
here (but see Reed et al. 2007, 2008). 

We perform this test here, using a large N-body simulation 
from the MICE collaboration (Fosalba et al. 2007) with 1024^ 
particles and cubic volume of side Lbox ~ 768 ' Mpc, for the 
cosmological parameters listed in §[T] The initial conditions 
were set at z = 50 using the Zel'dovich approximation, with an 
input linear power spectrum given by the analytic fit of Eisen- 
stein and Hu (1999). The subsequent gravitational evolution 
was followed using the Tree-SPH code Gadget-2 (Springel et 



al. 2005). Halos were identified using the Friends-of-Friends 
algorithm (Davis et al. 1985), with linking length equal to 
0. 164 times the mean interparticle density. The minimum halo 
mass resolved in the simulation is Mmin = 6 x 10"/?^' M©, 
with a minimum of 20 particles per halo. 

The halo mass function from the simulation is shown as 
solid symbols in Figure [T] (filled circles, filled squares, open 
triangles and open circles indicate the abundances of halos at 
redshifts z =3, 4, 5 and 6, respectively). We plot the quantity 
/(M,z) = nh{M,z)dM, where nh{M,z)dM is the number of 
halos per comoving volume at redshift z with mass between 
M and M + dM. The dashed and solid lines are the analytical 
models from Sheth & Tormen (1999; ST hereafter) and Jenk- 
ins et al. (2001; their equation B2), respectively, plotted at the 
same redshifts. To better display the difference between sim- 
ulations and models, the lower panel of Figure [T] shows the 
fractional deviation with respect to the Jenkins et al. (2001) 
fit, with all the lines and symbols as in the upper panel. 

Overall, we find that the ST model fits the simulations in the 
range z = 3 to 5 within 15% accuracy, while for z = 6 the error 
is about 20% in the range log(M/M0 = 12 - 13. We use 
the ST model in the rest of the paper, because the Jenkins et al. 
(2001) formula is clearly a worse fit to the simulation results 
in the regime of interest. 

As mentioned before, the halos we are interested in are rare, 
and studying their clustering properties is therefore difficult. 
This "rarity" can be quantified by means of the peak height 
V = Sc/(t{M,z), which characterizes the amplitude of density 
fluctuations from which a halo of mass M forms at a given 
redshift z (here, 6c ^ 1 .686, and (t(M,z) is the linear overden- 
sity variance in spheres enclosing a mean mass M). 

Gao et al. (2005) computed the halo bias using the Mille- 
nium simulation (Springel 2005) at redshifts z = — 5, but 
only for halos collapsing from fluctuations up to 3a. Angulo 
et al. (2008) measured the bias of ^ 4.5(t halos, but only for 
z < 3. In the context of the reionization of the universe. Reed 
et al (2008) studied the bias of < 4(t halos at redshift z > 10. 
Additional work on halo bias is presented in Seljak & Warren 
(2004), Cohn & White (2008), Basilakos et al. (2008), and 
references therein. In this section we extend these studies to 
the regime of our interest, namely 3 — 5 a halos at z = 3 — 6. 

We computed the bias factor of halos from simulation out- 
puts at z = 3,4,5,5.5 and 6. At each output we divided the 
halo catalogue into three mass bins of equal separation in 
logM, log(M/MQ/i-i) = 11.75- 12.25,12.25- 12.75, and 
12.75 — 13.25. We then measured the ratio of correlation 
functions /? = ^/„„(r)/^mm('") at 10 bins of equal width in logr 
in the range 8/!^ 'Mpc < r < 3 8/1^ 'Mpc (where ^/„„ is the 
two-point halo-matter correlation and f„„„ the matter-matter 
correlation function). The bias was computed as the mean of 
these values, and their variance was used as a rough estimate 
of the error. We warn that this error indicator may be underes- 
timating the true uncertainty in our measurements, since the 
correlation function errors are correlated in neighboring radial 
bins (although this effect is less severe in the presence of shot 
noise). 

We analyze halo bias from ^/„„ instead of to overcome 
the intrinsic noise in the latter quantity due to low halo abun- 
dance (see also Cohn & White 2008). The two definitions 
may differ owing to stochasticity in the halo-matter relation. 
However, we have tested that both measures yield consistent 
results (within error bars), while the variance among different 
bins is reduced by about 50% when using f/,„, (details of this 
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FlG. 2. — Halo bias estimated from ^i,„,/S,„im on scales 8 — 38/i~' Mpc~' . The symbols represent the results from a MICE simulation for different redshifts 
and halo masses as labeled (see text for details). In the left panel we show halo bias vs. mass and the coiresponding prediction of Sheth, Mo & Tormen (2001). 
The right panel shows bias vs. peak height u = 5c/cr and includes the Jing (1998) fit (which mostly coincides with Mo & White [1996] expression at these 
values of u). The Sheth et al. (2001) bias works well overall, but it underestimates the results from the simulation at high redshifts. Jing's fit, on the other hand, 
oveipredicts the measurements for all masses and redshifts studied by as much as 15 — 20%. 



comparison are given in the Appendix). 

In addition to the 8 — 38/j~'Mpc measurements, we have 
computed bias using the ^20 measure adopted by S07 and us- 
ing the range 30 — 60/;^' Mpc. These resuhs are reported in 
the Appendix. We find no evidence for scale dependence of 
the halo bias outside our statistical uncertainties, but the is- 
sue deserves further investigation in future work (Reed et al. 
2009). 

In Figure |2] we show the results for the halo bias at high 
redshift as obtained from the MICE simulation (with symbols 
corresponding to different redshifts as labeled in the figure). 
The left panel depicts the bias as a function of halo mass for 
various redshifts; the lines are predictions from the ellipsoidal 
collapse formula of Sheth, Mo & Tormen (2001), 

bsMT = IH 1=^ \y/a{av^) + ^/ab{av^y~'^ 
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(16) 



where a = Q.lQl, b = 0.5, c ^ 0.6, v ^ 5c/<t{M,z) and 5c = 
1 .686. The right panel shows instead the bias as a function of 
peak height v, in terms of which the predictions for all red- 
shifts coincide (equation |[T6l ). In addition to equation (fTST i, 
we also include in this figure the fitting formula derived by 
Jing (1998; see also Mo & White 1996). 

Figure |2]shows that Jing's (1998) fit overestimates the bias 
at all redshifts and masses studied at the 15 — 20% level. The 
Sheth, Mo, & Tormen (2001) prescription is in good agree- 
ment with the simulation for the lower fluctuations (j^ < 4) 
that correspond to halos of mass 3 x 10'^/;^' M© at z < 5, but 
it underestimates the bias of the rarest halos of > 4 by up 
to 10%. As noted by Crocce, Pueblas & Scoccimarro (2006), 
transients from the Zel'dovich dynamics generally used to set 
up the initial conditions lead to systematically high values for 
halo bias. This effect manifests itself more strongly in rare ha- 
los, so the discrepancy with Sheth et al. (2001) in this regime 



could be a numerical artifact rather than inaccuracy of the an- 
alytic model. Our conclusions are in good agreement with 
existing work on halo bias covering slightly different regimes 
(and mentioned at the beginning of this section). Therefore, 
we will use the Sheth et al. (2001) model for the bias and 
discuss in § 13.21 the impact that adopting Jing's (1998) bias 
formula would have on our conclusions. 

3. RESULTS 
3.1. COMPARISON TO OBSERVATIONAL DATA 
3.1.1. The data 

In this section we compare our model predictions with the 
available data on quasar clustering and the AGN luminosity 
function. This comparison is made in Figure 3, where the up- 
per left panel shows results for the AGN correlation length 
and the other three panels show the luminosity function at 
three different redshifts: z = 3.1, z = 4.5 and z = 6. 

The data on the clustering are taken from S07, who have 
recently extended beyond z ^ 3 previous measurements of 
the quasar clustering at lower redshifts from the Two Degree 
Field Quasar Redshift Survey (Porciani et al. 2004; Croom 
et al. 2005; Porciani and Norberg 2006; da Angela et al. 
2006; Mounthrichas et al. 2008) and SDSS (e.g., Myers et 
al. 2007; Strand et al. 2008; Padmanabhan et al. 2008). 
All the symbols in the upper left panel of Figure [3] represent 
the SDSS measurements averaged over sources with redshifts 
2.9 < z < 3.5 and z > 3.5. The diamonds and triangles refer to 
the S07 results extracted for the "good" and whole samples, 
respectively.^ Following S07, we take z = 3.1 and z = 4.0 as 

^ S07 in their clustering analysis remove the "bad" fields, i.e. those which 
do not fully satisfy their photometric criteria of completeness, but also re- 
port clustering measurements performed on the whole sample. We presume 
throughout this paper that the "good" measurements (shown with diamonds) 
are more rehable and the ones that any successful model must reproduce; 
however, following S07, we will always report both sets of data in the Fig- 
ures. The smaller number of pairs at small separations in the "good" sample 
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Fig. 3. — Upper left panel: Model prediction.s for the quasar correlation length ro as a function of redshift for different values of the input parameters, as 
labeled, computed above the luminosity threshold taken from Figure 1 of Richards et al. (2006a). The diamonds and triangles are the Shen et al. (2007) 
clustering measurements, corrected to Hq = 70km s~' Mpc~', calibrated on their "good" and total sample, respectively (see Shen et al. for details). Upper 
right panel: model predicted luminosity functions at z = 3.1, for the same set of models; the data are the collection from Shankar et al. (2009a) and Shankar & 
Mathur (2007), to which we refer the reader for details. Lower left panel: model predicted luminosity functions at z = 4.5. Lower right panel: model predicted 
luminosity functions at z = 6.0. The open and filled circles in the last three panels represent the AGN luminosity function before and after obscuration correction, 
respectively (see text for details). The vertical, thick, dotted lines in this and the following Figures mark the bolometric luminosity of L = 8 X 1 0"** erg s ^ ' , taken 
as the approximate luminosity threshold of the clustering measurements. Only data in the luminosity function above this threshold have been taken into account 
in the x^-fitting. 



the effective measurement redshifts for the two redshift bins. 

The S07 clustering measurements are for optically identi- 
fied AGNs only. However, the growth of black holes is con- 
nected to the total luminous output of the AGN population, 
not that of obscured or unobscured sources alone. We con- 
sider obscuration here as a random variable not linked with 
the large scale clustering of AGNs, so that the correlation 
length of obscured AGNs is the same as that of unobscured 
ones of the same bolometric luminosity. This assumption is 
plausible regardless of whether obscuration is principally a 
geometrical effect or an evolutionary phase. 

Following SWM, we take the AGN luminosity function at 
the mean redshifts of z = 3.1, 4.5, and 6, where most of the 
high redshift optical and X-ray data sets collected in SWM 
and Shankar & Mathur (2007; and references therein) are 
concentrated. We then adopt*" equation (4) of Hopkins et al. 
(2007a) to re-normalize the luminosity function to include ob- 
scured sources, assuming the obscuration is independent of 
redshift. However, because the obscuration correction may 

could lead to systematic errors in the correlation length estimate as well as 
larger statistical uncertainties (Y. Shen, private communication). 

' We also insert a Jacobian coirection factor in their equation (4) between 
observed B-band and bolometric luminosities. 



suffer from significant uncertainties, even up to a factor of a 
few (e.g., Ueda et al. 2003; La Franca et al. 2005; Tozzi et al. 
2006; Gilli et al. 2007; Hopkins et al. 2007a), when compar- 
ing model predictions to the bolometric luminosity function 
we will also include uncertainty in the obscuration correction 
as a source of systematic error to be added to the statistical 
error of the luminosity function measurements. In Figure [3] 
the filled and open circles show the luminosity function with 
and without obscuration corrections, respectively. 

3. 1 .2. General properties of the models 

Figures |3] and |4] illustrate the dependence of model predic- 
tions on the adopted parameters. In general terms, we can 
understand the interplay between the different parameters by 
combining equations (|5]i, (|9|l and ( fTTT i. Before examining the 
impact of individual parameter changes, we should note that 
there is one exact degeneracy within our family of models, 
if the luminosity function and correlation length are the only 
constraints. If we lower the Eddington ratio by a factor of F 
but raise the Mbh -Kir normalization a by the same factor, 
then the host halo mass at a given quasar luminosity is un- 
changed, so the predicted clustering is unchanged. All black 
hole masses are larger by a factor F, and so are their average 



7 



growth rates required to match the evolving halo mass func- 
tion, but if we lower the efficiency factor / = e/(l — e)by the 
same factor F, then the luminosity function implied by this 
growth is unchanged. Our analysis therefore cannot constrain 
A and / individually, but it can provide interesting constraints 
on the ratio fl A. 

In Figure [3] all models have evolution parameter 7=1.0 
and a tight correlation between Mbh and M, with E = 0. 1 dex. 
Solid lines in each panel show the predictions of a "reference 
model" with radiative efficiency e = 0. 15, Eddington ratio A = 
0.25, and a normalization of the Mbh -Kk relation a = 1.1. 
This model matches the SO? value of ro at z = 3.1. At z = 4, 
it is consistent with S07's measurement from the full quasar 
sample, but it falls below the "good" sample measurements by 
about 2(7. This model is in fairly good overall agreement with 
the bright end of the AGN luminosity function at all redshifts, 
though it is somewhat low at z = 4.5. 

In general all our models tend to overpredict the faint end of 
the AGN luminosity function below 0"*^ erg s ^ ' at z = 3 . 1 
and, more severely, at z = 4.5. These behaviors suggest that 
one or more of the model assumptions break down at lower 
luminosities. For example, the assumption of a constant A 
and € may not be valid. Alternatively, the assumed monotonic 
relation between black hole mass and halo mass could break 
down in this regime (see, e.g., Tanaka & Haiman 2008). How- 
ever, these hypotheses cannot be tested with the present data 
because the bias measurements by S07 do not probe luminosi- 
ties fainter than L < 10'*^ergs^', which is where our models 
start diverging from the data. We therefore do not attempt to 
reproduce the faint end of the AGN luminosity function with 
our models in this work. 

We now examine the consequences of varying each of the 
five model parameters, as listed at the end of § 12.11 We first 
consider lowering the Eddington ratio to A = 0. 1, keeping the 
other parameters fixed. Since the black hole abundances and 
their growth rates are fixed by their correspondence to ha- 
los, the duty cycles must increase as the inverse of A to com- 
pensate for the lower accretion rates during the active phase, 
thereby keeping the average volume emissivity from quasars 
constant. The results for this case are shown as the dotted line 
in Figure [3] A better match to the high observed clustering 
amplitude is clearly achieved, because the observed quasars 
correspond to more massive black holes and rarer halos. How- 
ever, the fit to the luminosity function is worse because the 
abundance of the most luminous quasars is underpredicted. 
Low values of the Eddington ratio are also disfavored by other 
observational (e.g., McLure & Dunlop 2004; Vestergaard et 
al. 2004; Bentz et al. 2006; Kollmeier et al. 2006; Kurk 
et al. 2007; Netzer & Trakhtenbrot 2007; Shen et al. 2008) 
and theoretical studies (Shankar et al. 2004; Lapi et al. 2006; 
Volonteri et al. 2006; Li et al. 2007; SWM; Di Matteo et al. 
2008). 

When the Mbh -Kir normalization a is lowered (dashed 
curve in Figure [3]l, the effect is simply to lower the black 
hole masses and quasar luminosities at fixed abundance. At 
fixed quasar luminosity, the clustering increases owing to 
the greater mass of the associated halos, but the decrease 
in abundance prevents a good match to the data. The dot- 
dashed curve illustrates the triple degeneracy described at the 
beginning of this section: a different set of (A,Q;,e) values 
whose predictions are nearly identical to those of the refer- 
ence model. 

Figure |4] shows the effect of varying the scatter S or the 
evolution parameter 7 of the Mbh -Kir relation. A low scat- 



ter maximizes the bias for a given set of other parameters, so 
low scatter is favored to reproduce the high values measured 
for ro at these redshifts (see also WMC). On the other hand, 
semi-empirical studies and AGN theoretical modeling support 
a significant intrinsic scatter for the Mbh -M relation at red- 
shifts z < 3 (e.g., Lapi et al. 2006; Haiman et al. 2007; Myers 
et al. 2007; Gultekin et al. 2009). Increasing the scatter to 
S = 0.3 (dotted lines), boosts the AGN luminosity function 
by increasing the number of massive black holes, but it de- 
presses clustering because more quasars at a given L reside in 
less massive halos. It also slightly flattens the dependence of 
the predicted ro on redshift. We then need to lower a or A 
to restore the luminosity function and increase ro. However, 
lowering A or a would also require higher radiative efficien- 
cies to keep the match to the luminosity function. We also find 
that models with E = 0.3, e > 0.25, and a < 0.7 predict more 
AGNs than black holes, yielding the unphysical condition of 
Po > 1 at z > 4. 

The parameter 7 regulates the amplitude of the model AGN 
luminosity function at high redshifts relative to that at z = 3 . 1 . 
Dashed lines in Figure |4] show a model with evolution index 
7 = and other parameters the same as those of the reference 
model. Lowering 7 maps the same L to higher mass, more 
biased halos at higher redshifts, steepening the ro — z relation 
and bringing it closer to the observed trend. However, more 
massive halos are rarer at higher redshifts, so the predicted 
AGN luminosity function drops significantly below the data 
at z = 4.5 and, especially at z = 6. Reproducing the observed 
luminosity evolution requires positive evolution (7 > 0) of the 
Mbh -Kir relation. 

3.1.3. A closer comparison 

Figure|5]presents a more systematic view of the dependence 
of clustering and luminosity function predictions on model 
parameters. Because none of our models reproduce all as- 
pects of the data, and because observational errors may in 
some cases be dominated by systematic rather than statistical 
uncertainties, we have taken only a semi-quantitative route to 
comparing models and measurements. The upper left panel 
shows models with A = 0.25, E = 0.1, and 7=1, defining 
the contour levels of acceptable models on a grid of {a, e) 
values. The blue and red areas define the regions where the 
Xjjjf for the luminosity is below 3 and 1.5, respectively. Here 
Xdof ~ X^/N, where = 45 is the number of points in the 
luminosity function, which include only those points (from 
Figure [3]) with L > 8 x lO'^^ergs^', the approximate lumi- 
nosity threshold of the clustering measurements, marked with 
vertical, thick, dotted lines in the Figures. For comparison, 
note that the models shown by the dashed and dotted lines 
in Figure[3]have x^/N = 5.41 and 18.13, respectively, while 
the reference model has x^/A^ = 1 .37. If the data points were 
independent, then even x^/N = 1.37 for = 45 would be 
an enormous statistical discrepancy, but the systematic uncer- 
tainty in the obscuration correction, at least, is highly corre- 
lated among points at a given redshift, motivating our rather 
loose criterion for "agreement". We assign observational er- 
rors to each data point equal to the reported statistical error 
(usually derived from the Poisson error on counts in the bin) 
summed in quadrature with 50% of the difference between 
the obscuration corrected and uncorrected luminosity func- 
tion estimates. This procedure is ad hoc, but it captures the 
reasonable expectation that the uncertainty on the obscuration 
correction is of the same order as (but smaller than) the correc- 
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Fig. 4. — Impact of changing the scatter E or evolution parameter 7 of the Mbh -K-ir relation. The format is as in Figure|3] and the adopted model parameters 
are labeled in the upper left panel. Increasing scatter worsens the match to the clustering data, and lowering 7 worsens the match to luminosity function evolution. 



tion itself, and the fact that the scatter among data sets visible 
in Figure[3]is comparable to the difference between open and 
filled symbols. 

The double-hatched and hatched areas define the regions 
where = ('"o.obs - '"o,pi-ed)^/o-obs, with the S07 value of 
(''o.obsjO'obs) — (24.3, 2. 4)/?"' at z = 4.0, is above 6 (i.e., a 
> 2.5(7 discrepancy) and 4, respectively. Note that the con- 
tour plots for the clustering are vertical, given that the pre- 
dicted clustering strength is independent of the values for the 
radiative efficiency (see §|2l). We find the constraints from the 
S07 clustering measurement of (ro.obs,crob,s) = (16.9, 1.7)/;^' 
at z = 3 . 1 not to be very constraining given that almost all 
models explored in Figure |5] are consistent with such data at 
the < 2(T level. It is clear from Figure |5] that the acceptable 
models are in general places in the upper-left corner of the 
(a, e) plane, characterized by higher radiative efficiencies and 
lower values of a for the same quasar luminosity /black hole 
mass, which implies higher halo masses (equation 15]) and 
corresponding clustering amplitude. 

Examination of Figure |5] reinforces the generality of the 
points made in our discussion of Figures |3] and |4] The cir- 
cle in the upper left panel marks our reference model with 
a = 1.1, 7 = 1.0, E = 0.1, A = 0.25, e = 0.15. Note that 
our reference model is not the best-fit model, as some mod- 
els characterized by higher radiative efficiency and lower a 
have an overall lower x^. However, we preferred to adopt 
as working models those defined by not too extreme values 
of the radiative efficiency. Also, the reference model already 



predicts /"o ~ 1 at z = 6 (see Figure [8] below), and lowering 
a reduces the black hole space density and pushes Pq above 
unity. Other models with the same a but different e have 
identical clustering, but the match to the observed luminos- 
ity function becomes worse for e < 0.1 and e > 0.25. Lower- 
ing a at fixed e improves the clustering agreement but quickly 
makes the luminosity function agreement worse. Raising a 
to 1.2 or 1.3 slightly improves the luminosity function agree- 
ment but worsens the clustering agreement. Raising A to 0.5 
(upper right panel) worsens the agreement with the z = 4 clus- 
tering if e and a are held fixed. However, because of the 
3-way degeneracy noted at the beginning of this section, a 
model with A ~ 0.5, e ~ 0.25, and a = 0.5 makes very sim- 
ilar predictions to a model with A = 0.25, e = 0.15, a = 1.0 
(which has //A smaller by a factor of « 2), and we disfavor 
the higher A models only on physical grounds because of the 
high required efficiency. Models with E = 0.3 (lower left) 
yield consistently worse agreement with the z = 4 correlation 
length unless lower values of a are adopted, but low-a mod- 
els produce Pq > 1 at high redshifts and require high radiative 
efficiencies to match the luminosity function. Models with 
7 = 0.0 (lower right) yield consistently worse agreement with 
the luminosity function. 

Our reference model underpredicts the S07 z = 4 corre- 
lation length (the value for "good" fields) by 2.2(7, and it 
slightly underpredicts the observed luminosity function in the 
luminosity range corresponding to the S07 quasar sample. If 
we take these discrepancies as a maximal allowed level of 
disagreement, then our reference model effectively defines a 
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Fig. 5. — The Xiof P^r degree of freedom as a function of the radiative efficiency e and a, the normahzation of the Mbh -Vvii relation, with other parameters 
fixed at the values listed on top of each panel. The blue and red areas define the regions in the e-a plane where the Xi„f the luminosity is below 3 and 
1.5, respectively. For the luminosity function we have used only the data with L > 8 X lO'^'ergs^', which is the luminosity threshold above which clustering 
measurements are available. The double-hatched and hatched areas define the regions where the Xj^f for the correlation length at z - 
respectively. The circle in the upper left panel marks the parameters of our reference model. 
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lower limit on the allowed value of //A, at //A=0.7. This 
conclusion does not depend on our adopted bolometric cor- 
rection. If we assumed a bolometric correction higher by a 
factor r, then we would require higher A (by the same factor) 
for fixed black hole masses to match our revised estimate of 
the bolometric luminosity function. We would also require 
higher /, again by a factor F, to reproduce the observed lu- 
minosity function history while building the same black hole 
population. For observationally estimated Eddington ratios 
A > 0.25 (Kollmeier et al. 2006; Shen et al. 2008) our limit 
on //A implies e > 0.15, significantly higher than the radia- 
tive efficiency e ~ 0.1 expected for the disk accretion onto a 
non-rotating black hole. 

3.1.4. Varying the bias 

These constraints would be much looser if we adopted the 
Jing (1998) bias function instead of the Sheth et al. (2001) 
formula that fits our N-body data. As akeady discussed in 
the previous sections, the Jing (1998) formula predicts a sig- 
nificantly higher value of the bias. Therefore, a much larger 
family of models can match the z — 4- S07 clustering measure- 
ments, with no strict requirement for a high //A ratio. For 
example. Figure |6] compares the predictions of the reference 
model to two alternatives, one with e = 0.065 (//A « 0.28) 
and one with A = 0.5 (//Aw 0.22), with ro calculated using 
the Jing (1998) formula in all models. The luminosity func- 
tion predictions of the reference model are unchanged, and 



all three models yield acceptable agreement with the z ~ 4- 
clustering measurement. The low e model underpredicts the 
luminosity function, but the A = 0.5 model overpredicts it, 
and lowering e to ^ 0.1 in this case would yield acceptable 
agreement. All three models overpredict the z = 3.1 correla- 
tion length. With optimal choices of a and A, one could find 
models that graze the top of the z = 3 . 1 error bar and the mean 
of the z = 4 correlation length while acceptably matching the 
bright end of the luminosity function (e.g., A = 0.3, e = 0.06, 
"f — 1, a = 2.3, S = 0.1). Alternatively, one could adopt any 
of the models shown in Figure |6] but drive down the z = 3.1 
clustering by assuming that the scatter E grows substantially 
between z = 4 and z = 3. 

3.2. BIAS AND DUTY CYCLE PREDICTIONS FOR THE 
REFERENCE MODEL 

Here we discuss further properties and predictions of a 
model that simultaneously matches the observed luminosity 
function and (at the 2a level) the clustering. For simplicity 
all the results presented below are obtained from the refer- 
ence model, which has (a, S, 7, A, e)=(l.l, 0.1, 1.0, 0.25, 
0.15). Figure I2] shows the predicted bias b{L,z) as a function 
of B-band magnitude; the solid, dashed and dotted lines corre- 
spond to z = 3.1, 4.5 and 6, respectively. Note that b{L,z) now 
refers to bias at a given luminosity rather than above a given 
luminosity. The predicted bias increases significantly with lu- 
minosity, at variance with what has been observed at lower 
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Fig. 6. — Same format as Figure[3] but Jing's (1998) bias formula has been adopted instead of the Sheth et al. (2001) one, and we consider models with lower 
values of e/X in addition to the reference model (solid line). In the upper left panel, the solid line overwrites the dot-dashed line because the two models have 
the same black hole mass-halo mass relation. 



redshifts z < 2, where evidence for a much flatter behavior of 
the bias against luminosity has been found (e.g., da Angela et 
al. 2006; Myers et al. 2007; Porciani & Norberg 2006; Coil et 
al. 2007; Mountrichas et al. 2008; Padmanabhan et al. 2008). 
Note that our prediction applies only to the very bright end of 
the AGN luminosity function, with L > 8 x 10"*^ erg s^' ; there 
are no available clustering measurements below this luminos- 
ity atz > 3.5. 
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Fig. 7. — Predicted bias for our reference model as a function of B-band 
magnitude at different redshifts, as labeled. 

We compute the effective halo mass that corresponds to 



the quasar hosts at the S07 luminosity thresholds via the re- 
lation b{Msff,z) = b{z)- Our reference model yields Mgff ~ 
1.1 X 10'^h"'Mo, nearly constant within 3 < z < 6. This 
mass scale is in marginal agreement with what has been in- 
ferred from the clustering analysis of large AGN-galaxy sur- 
veys at lower redshifts (e.g., Myers et al. 2007; Mountrichas 
et al. 2008), supporting a roughly constant host halo mass for 
luminous quasars at all times. Our reference model underpre- 
dicts the S07 correlation length at z = 4, so if we used their 
measured bias we would obtain a somewhat higher Meff at 
this redshift. In fact, S07 find an effective host halo mass for 
quasars at z > 3.5 that is a factor of two higher than the host 
halo mass for quasars with 2.9 < z < 3.5. Their mean values 
are Meff = 2.5 x IO'^/i-'Mq and5 x IO'^/j-'Mo in the low 
and high redshift bins, respectively, significantly lower than 
our quoted value, owing to their use of the Jing (1998) bias 
formula, which yields lower halo masses at fixed bias. 

Francke et al. (2008) have recently measured the clustering 
of 58 X-ray selected AGNs at z ^ 3 in the Extended Chan- 
dra Deep Field South, cross correlating them with a sample 
of 1385 luminous blue galaxies at the similar redshifts. Their 
quoted bias is b ~ 4.7 ±1.7 corresponding to halos of mass 

logM/M© — l2.6t_Q'l, i'^ li'^^ ^^^^ '■^^ previous findings by 
Adelberger & Steidel (2005) derived from optical quasar sam- 
ples at similar redshifts and luminosities. These studies probe 
AGNs about 6 magnitudes fainter than those probed by S07, 
and they seem to support a significant decrease of the bias at 
lower luminosities. These results would then be in qualitative 
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agreement with our model predictions, but at variance with 
the flat dependence of quasar clustering on luminosity found 
at lower redshifts (e.g., Myers et al. 2007; Coil et al. 2007; 
Padmanabhan et al. 2008). Larger surveys of low-luminosity 
quasars are needed to reduce the errors on the clustering mea- 
surements. 
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Fig. 8. — Predicted duty cycle for our reference model as a function of 
black hole mass at different redshifts, as labeled. The vertical dot-dashed line 
marks the point below which the results should be treated with caution, since 
there are no clustering constraints below this limit and the model overpredicts 
the luminosity function. 

Figure [8] shows the duty cycle Pq as a function of black 
hole mass and redshift (eq.[T2]). Results below ^ 6 x 10*^ M©, 
marked with a vertical dot-dashed line in the Figure, should be 
treated with caution, since there are no clustering constraints 
in this regime and the model overpredicts the luminosity func- 
tion. Above this limit the predicted duty cycles are roughly 
constant with mass, with values of Pq ~ 0.28,0.52 and 0.95 
at z = 3.1, 4.5 and 6, respectively. We are using the model 
luminosity function to compute these duty cycles via equa- 
tion ( fT2] ). but the agreement with the observed ^{L,z) is good 
enough (Figure O that we can consider this a smooth proxy 
for the observational data. 

4. PREDICTIONS FOR Z > 6 

The high duty cycle inferred at z = 6 has profound implica- 
tions for the evolution of the luminosity function at still higher 
redshifts. Between z — 'i and z = 6, the decreasing abundance 
of halos with increasing redshift is partly compensated by the 
factor of three increase in duty cycle. However, duty cycles 
cannot exceed unity by definition, so at z > 6 the fast drop 
of the massive and rare host halos implies an equally rapid 
decline in the number density of luminous quasars. At the 
same time, the implied mass of quasar hosts moves even fur- 
ther out on the exponential tail of the halo mass function. Our 
models thus predict a decline in high redshift quasar numbers 
much steeper than expected from simple extrapolations of the 
z = 3 — 6 luminosity function. 

Figure |9] demonstrates this point, showing the reference 
model predicted number counts of AGNs per square degree 
per unit redshift as a function of redshift, above luminos- 
ity thresholds of logL/ergs^' =47,47.5, and 48, as labeled. 
Evolution is more rapid for higher luminosity AGN because 
their host halos are further out on the tail of the mass func- 
tion. The thin lines refer to extrapolations of the Fan et al. 
(2004) luminosity function at the same redshifts and lumi- 
nosities. The latter is a power law <i>(L) oc L^^ ' that describes 
the statistics of optical quasars in the range logL/erg s^' ^47 
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Fig. 9. — Number counts of AGNs per square degree per unit redshift as 
predicted by our reference model (thick lines) as a function of redshift, above 
the labeled luminosity thresholds. Thin lines refer to extrapolations of the 
Fan et al. (2004) luminosity function at the same redshifts and luminosities. 
The large open squares indicate the number density per unit redshift con'e- 
sponding to one single observable quasar in the whole sky. 



and 5.5 < z < 6.5; we do not apply any obscuration correction. 
Fan et al. (2004) find that a good representation of the data 
requires redshift evolution $(L,z) oc 10^""^**~. It is evident 
from Figure |9] that our reference model predicts a decrease in 
AGN number density much faster than the one expected by 
naively extrapolating the Fan et al. (2004) trend to z > 6. The 
open squares in Figure |9] indicate the number density per unit 
redshift corresponding to one single observable quasar in the 
whole sky. According to this model, the highest-z quasar in 
the sky with true L > lO"*^'^ erg s^' should be at z ^ 7.5 in our 
model; all quasars detected with higher apparent luminosities 
by future surveys would have to be magnified by lensing (see 
also Richards et al. 2004). 

As recently discussed by Fontanot et al. (2007), the sur- 
face density inferred from the luminosity function of Fan et al. 
(2004) and Shankar & Mathur (2007) predicts that only a few 
luminous sources will be detected in the field of view of even 
the largest and deepest future surveys such as JWST, and EU- 
CLID. Our predictions suggest that these detections will be 
even rarer than simple empirical extrapolations predict. While 
the predictions of Figure |9] are specific to our adopted model 
parameters, this conclusion is likely to apply more generally 
to models that reproduce the strong z = 4 clustering found by 
S07. This clustering implies high host halo masses and hence 
high duty cycles at z = 4, so the declining black hole mass 
function cannot continue to be compensated by higher duty 
cycles towards higher redshifts (though rapid evolution of the 
Mbh -Vvir relation or rapidly increasing A values could com- 
pensate in principle). Very similar results are found with the 
e = 0.1 model of Figure[3] 

Because the host halos of high redshift quasars are so 
highly biased, the predicted clustering remains strong at z > 
6. For the e = 0.15 reference model, the predicted correla- 
tion length ro as a function of B-band luminosity and red- 
shift can be well approximated by the relation r()(MB,z) — 
27 X [(1 +z)/7]'"(-26.5+Mb)°'5 Mpc, valid in the range 
-29 <Mb< -26.5 and 6 < z < 9. 

Local observations imply a ratio Mbh /A^star ~ 1 .6 x lO^-' 
between the mass of the black hole and the stellar mass of its 
host bulge (e.g., Mai'coni & Hunt 2003, Haiing & Rix 2004). 
As discussed above, our reference model predicts increasing 
black hole masses at fixed virial velocity (7 > 1) at z > 3 as re- 
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quired to match the number density of very luminous quasars 
of z = 6. However, with the assumption that the baryon frac- 
tion within a halo is universal, this implies that an increas- 
ing fraction of the baryons must be locked up into the central 
black hole. We show this in Figure [TOl which plots the black 
hole-to-baryon fraction as a function of redshift predicted by 
our reference model. The baryonic mass fraction within any 
halo is set to be Mbar/^ — ft — 0.17 (e.g., Spergel et al. 
2007; Grain et al. 2007). Given that the stellar mass Mstar in 
the local universe is unlikely to exceed fh, and is more typi- 
cally < ft/3 in massive galaxies (e. g., Shankar et al. 2006), 
this implies that the black hole mass is < 1.6 x 10^^/3 the 
mass of the total baryons in the host halo. This local ratio be- 
tween black hole and baryonic mass is shown as a horizontal 
dotted line in Figure [TO] The fraction of baryons locked in 
the central black hole increases at higher redshifts following 
the increase of virial velocity at fixed halo mass (Eq. ID) and 
the increase of black hole mass at fixed virial velocity propor- 
tional to (1 (Eq Is)). For our reference model, the Mbh 
/Mbar ratio grows rapidly at high redshifts and exceeds the 
local value by nearly an order of magnitude at z > 6 for Mbh 
> IO^Mq. Note that even a model with 7 = still produces 
an increase of the Mbar/M ratio with redshift, driven by the 
redshift dependence in equation (HJl, although it is just a factor 
of a few in this case. 

0.1000 



0.0100 



0.0010 




0.0001 



Fig. 10. — Ratio between black hole and baryon mass within the halo, the 
latter computed as Mbar = 0.17 x M, for three values of the black hole mass, 
as labeled. This ratio at z > 4 gets higher than the local value between black 
hole and bulge mass of (1/3) X 1.6 X 10~' (dotted line; see text), implying 
that at fixed stellar mass, a larger fraction of the baryons in high mass halos 
is locked in the central black hole at early times. 



Therefore, we conclude that the relation between black hole 
and spheroidal stellar mass determined locally cannot con- 
tinue to hold at very high redshifts if the large clustering 
strength reported at z = 4 is to be matched, and that a much 
larger fraction of baryons in galaxies must accrete to the nu- 
clear black holes at z > 4. 

5. GOMPARISON TO GONSTANT DUTY GYGLE 
MODELS 

The results in § |3] show that matching the high clustering 
signal measured by S07 requires a high duty cycle Pq, which 
corresponds to quasars preferentially residing in high mass, 
less abundant halos. This result has also been discussed by 
S07 and by WMG, following the method outlined by Martini 
& Weinberg (2001) and Haiman & Hui (2001). The model de- 
scribed in §|2]assumes an a priori relation between luminosity 
L and halo mass M. Since this model also predicts the AGN 



luminosity function from the equation governing the growth 
of black holes, it implicitly predicts the duty cycle required 
to assign an AGN luminosity to a halo mass and match their 
abundances. Martini & Weinberg (2001) and WMG instead 
define the relation between L and M a posteriori, i.e., from the 
cumulative matching between the observed AGN luminosity 
function and the halo mass function, once an input duty cy- 
cle has been specified. Since both methods assume a (nearly) 
monotonic relation between luminosity and halo mass, they 
should yield a similar connection of duty cycle and cluster- 
ing in cases where the a priori model matches the observed 
luminosity function. 

To compare the two approaches in detail, we compute the 
relation between black hole mass and virial velocity for fixed 
duty cycle via the equation 

$BH(>MBH,z) = ^^^^=$/,(>Vvir,z), (17) 

where $(> L,z) is the model predicted AGN luminosity func- 
tion, and ^it{Vwii,z) is derived from the halo mass function 
and the Vvt -M relation of equation (|4|i. We assume A — 0.5 to 
convert from Mbh to L. Figure[TT]plots the relations implied 
by equation ( fTTI ) at redshifts z = 3.1 and z = 4. Gurves from 
top to bottom show the Mbh -Kir relation assuming a constant 
duty cycle Pq — 0.1,0.2,0.3 and 0.5. Higher Pq corresponds 
to rarer halos, hence higher Vyt and stronger clustering. The 
solid curves in the two panels represent the output Mbh -Kir 
relation for our reference model, which predicts a duty cy- 
cle of Po(z = 3.1) - 0.30 and Poiz = 4.0) - 0.50 at the high 
mass end. As expected, the Mbh -Kir relation of our reference 
model is similar at high masses to that model with similar duty 
cycle. At lower masses, our model does not perfectly match 
the observed luminosity function and does not predict a con- 
stant duty cycle. 

We compute the average bias for the constant duty cycle 
model via 



b{z) = 



f2^^^^^b{M,z)<f„iM,z)d\ogM 



iZ.iz) ^h{M,z)d\ogM 



(18) 



where Mminfz) is the halo mass corresponding to Lmin (z) via 
equations ^ and ( [TtT i. Because the relation between L and 
K'ir is determined by matching space densities, the predicted 
bias is independent of A. Figure [12] plots the corresponding 
correlation length ro, computed through equation ( fT3] ), as a 
function of Pq at z = 3. 1 and 4.0. Solid and dashed lines show 
the results of using the Sheth et al. (2001) and Jing (1998) 
formulas, respectively, for the bias b{M,z) in equation (fTSb . 
Shaded regions show the la range of the S07 measurements. 
As expected the clustering strength strongly increases with 
increasing duty cycle. In agreement with WMG, we find that 
a duty cycle Pq > 0.2 is required to reproduce S07's z = 4 
measurement with the Jing (1998) bias formula. However, our 
N-body results favor the Sheth et al. (2001) bias formula, and 
in this case we cannot match the S07 "good" measurement 
within 1 (T at z = 4 even for a maximal duty cycle Pq— I. 

As shown in § [3] models with high duty cycles require 
high //A ratios to reproduce the observed luminosity func- 
tion. This connection can be simply understood from equa- 
tion ( fTTI ): increasing the duty cycle decreases the number den- 
sity of halos that host AGNs, which in turn need to increase 
their e-folding time fgf to maintain the same observed lumi- 
nosity density. Figure [13] shows this effect in detail. We first 
integrate our reference model luminosity function from z 6 
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Fig. 1 1 . — Correlation between black hole mass and halo virial velocity implied by the cumulative matching between the halo and black hole mass functions, 
the latter derived from the reference model luminosity function and a constant input duty cycle Pq. The different Unes refer to different values of the duty cycle, 
as labeled, while the solid line is the Mbh -Vy\i relation con'esponding to our reference model. The left and right panels show the resulting relations at redshifts 
z = 3.1 and z = 4, respectively. 



down to a given redshift z as 



Pbh(> logL,z) 



1 



logL 



dt 



dz' 



dlogL' . 
(19) 

We consider only luminous AGNs that shine with luminosity 
logL/erg s^' > 45, corresponding to black hole masses above 
Mbh ^ IO^Mq, which ensures that we are properly tracking 
the accretion histories of the most massive black holes. It 
is evident from equation (T% that the accreted mass density 
does not depend on the black hole Eddington ratio distribution 
but only on the radiative efficiency. Our results are shown as 
stripes in Figure [13] for three different values of the radiative 
efficiency e = 0.10,0.15 and 0.25, from top to bottom. The 
left and right panels of Figure [13] show the integrated mass 
density at z = 3.1 and 4, respectively. Note that these are 
the black hole mass densities implied by the Soltan (1982) 
argument given an input luminosity function that is a good 
match to observations. 

Alternatively, by assuming an average duty cycle Pq{z) at a 
given redshift z we can convert the AGN luminosity function 
into a black hole mass density via 

f°° L' 

p'^^{>\ogL,z)^ —<f{L',z)dlogL'. (20) 
JiogL 

The latter estimate^ depends inversely on the Eddington ratio 
A because of the mapping of L to Mbh , but it does not depend 
on the radiative efficiency. We plot pgjj as a function of the 
duty cycle Pq for A = 0.25, 0.5, and 1.0, with solid, dashed, 
and dotted lines, respectively. Results for the z = 3 . 1 and z — 4- 
accreted mass density are shown in the left and right panels of 
Figure [13] respectively. It is noteworthy that the high radia- 
tive efficiency of e = 0.15, as used in our reference model, is 
consistent with Po{z = 3.1) 0.28 and Po{z = 4) 0.50, in 
perfect agreement with our findings presented in Figure [8] 

Overall, we find evidence for a general rule of thumb: if 
black holes accrete at a significant fraction of the Eddington 
luminosity (A > 0.25) and possess high duty cycles as derived 
from their strong clustering (Figure [T2T i. then they must also 

' Note that we do not cons ider any scatter between black hole mass and 
AGN luminosity in equation )20t , as the lumi nos ity function has been de- 
rived from the continuity equation in equation which requires a strictly 
monotonic relation between black hole and halo mass. 



radiate at high radiative efficiencies (e > 0.15) to match the 
AGN luminosity function and its evolution with redshift. This 
conclusion from constant duty cycle models is entirely consis- 
tent with our conclusions from Mbh -Kir models discussed in 
§§|2]|4] 

6. DISCUSSION AND CONCLUSIONS 

We have investigated constraints on the host halos, radia- 
tive efficiencies and active duty cycles of high redshift black 
holes that are implied by recent measurements of the AGN 
luminosity function at 3 < z < 6 and of optical quasar clus- 
tering at z « 3 and z ~ 4. In this work we have derived the 
predicted AGN luminosity function implied by a model black 
hole mass function. The latter is built from the dark matter 
halo mass function at each redshift by applying a model re- 
lation between black hole mass and halo virial velocity, mo- 
tivated by local observations. Our models are parameterized 
by the high-redshift normalization a and redshift evolution 
index 7 of the mean Mbh -Kir relation (equation ||3]), by the 
log-normal scatter S about this relation (in dex), and by the 
Eddington ratio A and radiative efficiency e of black hole ac- 
cretion. 

A reference model with (a,7,I],A,e) — 
(1.1,1.0,0.1,0.25,0.15) provides a good fit to the z = 3 
correlation length ro and a reasonable fit to the bright end 
of the luminosity function (L > 10'*^'^ergs^') at z = 3 — 6. 
It overpredicts the faint end of the luminosity function, 
probably indicating that our assumption of a constant A or 
power-law Mbh -Kir relation breaks down in this regime. 
More significantly, the model prediction is below S07's 
estimate of ro for luminous quasars at z = 4, by about 2a. 
While lowering a or A raises the predicted ro, it lowers 
the predicted luminosity function below the observations, 
unless we allow efficiencies greater than e = 0.15. Increasing 
the scatter S reduces the predicted clustering, making the 
overall fit to the data worse. If we use S07's "all field" 
estimate of ro instead of their "good field" estimate, then 
the discrepancy at z = 4 is under la. The reference model 
predicts substantial luminosity and redshift dependence 
of the quasar correlation length at z > 3 (Figure O, with 
ro w 27 X [(1 4-z)/7]"-3(-26.5 +Mb)°-5 Mpc for 6 < z < 9. 

Models that successfully match the high redshift bias at z > 
3 require luminous AGNs to reside in massive and therefore 
rare halos, implying high duty cycles, Pq ^ 0.2, 0.5, 0.9 at 
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Fig. 12. — Predicted clustering correlation length rg computed above the minimum survey sensitivity as a function of duty cycle and adopting the luminosity 
function deiived from our reference model. The solid and long-dashed lines refer to the imphed by using Jing's (1998) and Sheth et al.'s (2001) bias formula, 
respectively. The left and right panels show our results at redshifts z = 3.1 and z = 4, respectively. Dark and Hght shaded bands show the Icr range of the S07 
measurements at these redshifts, from "good fields" and "all fields", respectively. Maximal values of the duty cycle predict a clustering strength only marginally 
consistent with the data at z = 4. 



z = 3.1, 4.5, 6.0 in our reference model. Note that, although 
this model is consistent with the z — 4- clustering only at the 2a 
level, it already produces a duty cycle close to unity at high 
redshifts. Raising the predicted correlation requires putting 
quasars in more massive, less numerous halos, and thus tends 
to push the required duty cycle above unity. 

To simultaneously reproduce the observed luminosity func- 
tion and bias, models must have //A > 0.7, where / = 
e/(l — e), so that the mass density of black holes in these rare 
halos corresponds to the cumulative emissivity of the lumi- 
nous AGN. These findings are robust against uncertainties in 
the obscured fraction of AGNs or in the precise value of the 
mean bolometric correction (see discussion in § 13.1.31 ). The 
underlying physics that leads to these findings is easy to un- 
derstand. The strong observed clustering at z = 4 implies a 
high duty cycle and thus a low space density of massive black 
holes. Reproducing the observed AGN emissivity with the 
low total mass density in black holes requires a high radiative 
efficiency. Lowering the assumed Eddington ratio implies a 
higher mass density (because each black hole is more mas- 
sive) and a proportionally lower /. As shown in the Appen- 
dices, mergers are expected to have little impact on the BH 
mass function at these redshifts, but to the extent they do have 
an impact they raise the limit on // A by adding mass without 
associated luminosity. 

For any choice of the mean Eddington ratio, our success- 
ful models require positive evolution of the Mbh -Kir relation 
(7 > 0) at z > 3 to reproduce the evolving bright end of the 
luminosity function. Evolution of the Eddington ratio itself 
(higher A at higher z) could in principle yield similar evolu- 
tion. 

It is beyond the purposes of the current work to extrapo- 
late the simple model outlined here to lower redshifts. First 
of all, the basic treatment presented by SWM has shown that 
the large scale clustering of quasars can be simply matched 
by accretion models which evolve the black hole mass func- 
tion assuming reasonable values of the radiative efficiency and 
Eddington ratios, which satisfy Soltan (1982)'s constraint. 
Moreover, at lower redshifts, several additional physical in- 
puts need to be added to the model (e.g., the fraction of ac- 
tive satellites, mass-dependent Eddington ratios, AGN feed- 
back) to reproduce the full quasar clustering at all luminosi- 
ties, scales, and redshifts (e.g., Wyithe & Loeb 2003; Scanna- 



pieco & Ho 2004; Hopkins et al. 2008; Bonoli et al. 2009a; 
Thacker et al. 2009). 

Previous work that attempted to simultaneously match the 
quasar luminosity function and bias has yielded somewhat 
different results from ours. Wyithe & Loeb (2003, 2005; see 
also Rhook & Haehnelt 2006) developed a model aimed at 
reproducing both the bias and the AGN luminosity function 
at several redshifts. They expressed the relation between the 
luminosity and halo mass via some AGN feedback-motivated 
models for the black hole-halo relation, and they assumed that 
black holes grow at the Eddington limit and radiative effi- 
ciency of e = 0. 1 . Their values of // A would then be lower by 
a factor of several with respect to ours. These differences are 
due to a different AGN bolometric luminosity function used 
(ours being a factor of a few higher) and the absence of the 
SDSS bias measurements at z > 3. In brief, we do not think 
that these models would reproduce the observational data con- 
sidered in this paper. 

Our lower limit on //A translates to a lower limit on ra- 
diative efficiency e > 0.7A/ (1 + 0.7A). With observationally 
estimated values A w 0.3 for the Eddington ratios of lumi- 
nous high-redshift quasars (Kollmeier et al. 2006; Shen et 
al. 2008), this limit implies e > 0.17. Using a different ap- 
proach that links the observed AGN luminosity function to 
the local black hole mass function via the continuity equa- 
tion, a differential generalization of Soltan's (1982) cumula- 
tive emissivity argument, SWM estimate an average radiative 
efficiency e w 0.05 — 0.10. Other authors pursuing similar 
approaches and adopting similar bolometric corrections have 
reached similar conclusions (e.g., Salucci et al. 1999; Mar- 
coni et al. 2004; Shankar et al. 2004; Hopkins et al. 2007a; 
Yu & Lu 2008; see SWM for further discussion). As dis- 
cussed in detail by SWM, uncertainties in the local black hole 
mass function, bolometric corrections, and obscured fractions 
still leave significant range in the inferred value of e, but these 
uncertainties would have to be pushed to their extremes to ac- 
commodate e > 0.17. 

One possible resolution of this tension is that the typical 
radiative efficiency is higher at z > 3 than it is at the lower 
redshifts that dominate the overall growth of the black hole 
population. However, the similarity of quasar spectral en- 
ergy distributions at low and high redshifts (e.g., Richards 
et al. 2006b) argues against a systematic change in accre- 
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tion physics. We should therefore examine the loopholes in 
the argument for high efficiency presented here, noting that it 
is above all the z = 4 clustering measurement from S07 that 
drives our models to massive, rare halos and thus to high ef- 
ficiencies to reproduce the luminosity function. Adopting the 
Jing (1998) bias formula instead of the Sheth et al. (2001) 
formula would allow us to match the clustering with less mas- 
sive halos, but our numerical simulation results show that the 
Sheth et al. (2001) formula is more accurate in the relevant 
range of halo mass and redshift. 

Our conclusion is insensitive to the specific assumption of 
a one-to-one power-law relation between black hole mass and 
halo virial velocity: monotonically matching luminosity func- 
tions and halo mass functions leads to a similar conclusion 
(^ WMC), and adding scatter, while plausible on physi- 
cal grounds, only reduces clustering and thus exacerbates the 
underlying tension. Because the implied characteristic halo 
mass is already well above M*, halos hosting two or more 
quasars should be far too rare to significantly alter the large 
scale bias. Small-scale clustering studies of large z > 3 quasar 
samples (Hennawi et al. 2009; Shen et al. 2009), will set more 
definite constraints on the actual fraction of active subhalos. 

Our modeling does assume that halos hosting quasars have 
the same bias as average halos of the same mass, while Wech- 
sler et al. (2006) find that the youngest 25% of high redshift, 
high M jM-j, halos have a correlation length ^13% higher than 
the average correlation length of halos at the same mass and 
redshift. Clustering could be slightly boosted if active quasars 
preferentially occupy these younger halos, but the high duty 
cycle required in our models effectively closes this loophole, 
implying that the quasar host population includes the major- 
ity of massive halos rather than a small subset. Wyithe & 
Loeb (2009) suggest that the strong z = 4 clustering could 
be explained by tying quasar activity to recently merged ha- 
los, which might have stronger bias (see also Furlanetto & 
Kamionkowski 2006; Wetzel, Cohn, & White 2009). How- 
ever, we suspect that halos with substantial excess bias might 
be too rare to satisfy duty cycle constraints, and a recent study 
by Bonoli et al. (2009b) uses outputs from the Millennium 
Simulation (Springel et al. 2005) to show that no excess bias 
is present in recently merged massive halos. 

Perhaps the most plausible loophole in our conclusion is 
simply that the S07 z = 4 correlation length is overestimated 
(or its statistical error underestimated), since it is the first mea- 



surement of its kind and there is a significant difference be- 
tween the S07 values for "all fields" and "good fields" (even 
though only 15% of fields are excluded from the latter sam- 
ple). Even the central value of the (lower) "all fields" mea- 
surement requires high e or low A when combined with lu- 
minosity function constraints, but the lower Icr value can be 
reconciled with e w 0.1 and A « 0.25. The DR7 SDSS quasar 
sample should afford substantially better statistics than the 
DR5 sample analyzed by S07, allowing stronger conclusions 
about the host halo population. 

In these models, the rapid decline in the number of lumi- 
nous quasars between z = 3 and z = 6 is driven by the rapidly 
declining abundance of halos massive enough to host them. 
However, the drop in halo abundance is partly compensated 
by a rise in the duty cycle over this interval, from ^ 0.2 to 
^ 0.9 in our reference model. Since duty cycles cannot exceed 
one by definition, this compensation cannot continue much 
beyond z = 6, and the decline in host halo abundance accel- 
erates because these halos are far out on the exponential tail 
of the mass function. The predicted evolution of the quasar 
population at z > 6 is therefore much more rapid than simple 
extrapolations of the observed z = 3 — 6 behavior. This break 
to more rapid evolution at z > 6 should be a generic predic- 
tion of models that reproduce the strong observed clustering 
at z = 4, though in principle it could be softened by a rapid 
increase of Eddington ratios at z > 6 or by a sudden change in 
evolution of the Mbh -Kir relation. Surveys from the next gen- 
eration of wide-field infrared instruments will have to probe 
to low luminosities to reveal the population of growing black 
holes at z > 7. 
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APPENDIX 
HALO BIAS AT HIGH REDSHIFT 

In this Appendix we provide additional details of our bias analysis from the N-body simulation introduced in § 12.21 in particular 
the difference between using quasi-linear scales (e.g. 8 — 38 Z;^' Mpc as reported) or larger ones. We also comment on the 
distinction between deriving the bias from the halo autocorrelation ^/,;, or the halo-matter correlation function ^/„„. We finish by 
comparing our results with those of S07. 

The simulation used was provided by the MICE collaboration (Fosalba et al. 2008) and contains 1024^ particles in 
a cubic volume of side Lbox ~ 768/z^'Mpc, with cosmological parameters 51,„ — 0.25, JIl = 0.75, g% = 0.8, n ~ 0.95, 
h = //o/lOOkms^'Mpc^' = 0.7 and fti, = 0.044. Halos were identified in the comoving outputs at z = 3,4,5,5.5 and 6 us- 
ing the friends-of-friends algorithm with linking length b — 0.164. Finally, at each redshift the corresponding halo catalogue 
was divided in three non-overlapping sub-samples ofhalo masses in the ranges [5.6—17.7] x 1O"M0, [1.7 — 5.6] x lO'^M© and 
[5.6 — 17.7] X IO'^Mq respectively (i.e. bins of equal width in logM). 

We have obtained the halo bias from the ratio of correlation functions £,hm{r) / £,mm{r) averaged over 10 bins of equal length 
in logr in the radial range 8/i^'Mpc < r < 38/i^'Mpc. We also implemented the same measurements in the radial range 
28/i^'Mpc <r< 60/z^'Mpc in order to test for any dependence of the bias on scale. In Fig. lAll we show the ratio ^hm/£,mm at 
both ranges for all redshifts and mass bins studied. On the one hand, this figure indicates that within the intrinsic scatter there is 
no significant scale dependence of the bias at smaller separations. On the other hand, the values of the measured bias rise by only 
2 — 3 % when using 8 — 38 /i^' Mpc instead of 28 — 60 /i^' Mpc, but this difference is within the variance of the simulation. 
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Fig. 13. — Left panel: the horizontal stripes show the integrated black hole mass density above L = 10*^ ergs^ ' and 3.1 < z < 6 derived from our reference 
model luminosity function for three different values of the radiative efficiency, as labeled; triple dot- dashed, solid, long-dashed, and dot-dashed lines indicate the 
black hole mass density implied by the luminosity function at z = 3.1 integrated over mass assuming a mean Eddington ratio A = 0.1,0.25,0.5, 1, respectively, as 
a function of duty cycle Pq. Right panel: same as left panel but integrating the luminosity function over the redshift range 4 < z < 6. For acceptable combinations 
of A, e, and Pg, the line corresponding to A should intersect the band corresponding to e at duty cycle Pq. High duty cycles require high radiative efficiencies or 
low Eddington ratios to reconcile the cumulative AGN emissivity with the black hole mass density in rare halos (see text). 



In order to overcome the shot noise due to the low abundance of rare halos we decided to measure the bias from cross corre- 
lating the halo distribution to that of the dark matter This allowed us to extend the measurements to cases where the halo-halo 
correlation is too noisy to define a meaningful bias. In Table 1 we report the bias results for \/S,hh/Cmn (left table) and for £,hml£.mm 
(right table). In both cases we measured at scales in the range 8 /i^ ' Mpc < r < 38 Z;^ ' Mpc. The reported error corresponds to 
the variance among different bins (which might be taken as a rough representation of the true error of the measurement with the 
caveat that different bins are correlated). We find consistent results for the values of the bias derived from the two methods in 
those bins of mass and redshift where a reliable estimate can be obtained. 



log(M/MQ) = 


11.75-12.25 


12.25-12.75 


12.75-13.25 




log(M/Mo) = 


11.75-12.25 


12.25-12.75 


12.75-13.25 


z = 3 


/) = 3.95 ±0.07 


fo = 5.3±0.17 


fo = 7.8 ±0.69 




z = 3 


fo = 3.97 ±0.06 


fo = 5.27 ±0.09 


fo = 7.38 ±0.2 


z=4 


/) = 5.97 ±0.22 


fo = 8.1 ±0.57 






z=4 


fo = 5.9±0.13 


fo = 7.88 ±0.22 


b= 11.52±0.54 


z = 5 


/:, = 8.44 ±0.61 


h= 11.6±2.78 






z = 5 


fo = 8.25 ±0.27 


h= 11.44 ±0.5 




z=5.5 


h= 10.1 ±1.1 








z=5.5 


fo = 9.53 ±0.39 


h= 12.78 ±1.4 




z = 6 


h = 12.3±1.6 








z = 6 


b= 10.96 ±0.69 







TABLE Al 

Halo bias obtained from h = ^/g;,/, /gmm {left TABLE) ORh = £,h,„/^mm {right TABLE). THE VALUES ARE CONSISTENT WITH EACH OTHER WITHIN 

THE BIN-TO-BIN SCATTER, WHICH IS REPORTED AS THE CORRESPONDING ERROR. 



Finally, we compare our results to those of S07, who obtained the halo bias from b = ^ ^20/^20" with, 

ii^^^ r'\{ry^dr, (Al) 

where rmin — 5/i^' Mpc and ^max = 20/i~'Mpc. Using all halos more massive than Mmin = 2 x IQ^'h^^ Mq, S07 finds bsfs{z = 
3) = 6.2 and b^if{z = 4) = 10.2 (respectively 17% and 5% lower than Jing's [1998] prediction). Using a similar mass-cut to 
S07 and measuring the bias in the same way (i.e. from ^20) we obtain b^fi{z = 3) = 6.07 and b^fi{z = 4) = 9.35 (where we have 
included a 6% correction due to the difference in the assumed cosmology). These values are in good agreement (within 8%). 
However, and contrary to S07, we have chosen Sheth, Mo & Tormen (2001) expression as our primary bias prediction for reasons 
already outlined in § 12.21 

THE CONTRIBUTION OF BH MERGERS TO MASS ACCRETION 

In this Appendix we compute the expected contribution from mergers to the overall mass growth of the central black hole of a 
halo with mass Mq '-^ 1 — 2 x 1O'^M0 at z ~ 3 — 4, typical of the z > 3 quasar hosts studied in this paper. 

We trace the average mass growth history M(z) of a such halos at any redshift z by imposing the number density conservation 
within a comoving volume 

n[M(z),z]=n[Mo,z = 3]. (Bl) 

The result is shown with a long dashed line in the left panel of Figure IB II which shows that such halos grow, on average, by 
about a factor of ^ 6 within the redshift range 3 < z < 6. 
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Fig. A1. — Halo bias as a function of scale. Right panels: the ratio §;„„/5„i„i at large scales for different halo masses and redshifts, as labeled. Left panels: 
the same ratio at smaller separations, that encompass the scales where AGN clustering has being measured by S07. At smaller separations there is not significant 
scale dependence while the scatter of the measurement is lower than for large separations. The measured bias is higher by 2 — 3 % in the left hand panels, but this 
difference is within the scatter 

We then compute the expected growth of such halos due to halo and subhalo mergers. To such purpose we estimate the average 
number of mergers per unit redshift dN/dz experienced by a given halo of mass M{z) at redshift z by integrating the halo merger 
rates 

dz n[M[z),z] 

with ^min =0.1 and f^ax = 1- We take the halo merger rates per unit halo B[M{z),£,,z]/n[M{z),z] from Fakhouri & Ma (2008), 
given as empirical fits to the Millennium Simulation, with n[M{z),z] the number density of halos of mass M{z) at redshift z- The 
quantity B[M{z),£,,z] is the instantaneous merger rate at redshift z of halos with mass M{z) in units of /i"* Mpc^^dz^' d^^'Mg', 
with ^ = M1/M2 the mass ratio between the masses of the progenitor merging halos with total mass M{z) — My +M2. The total 
mass accreted AM (z) on the halo of mass M{z) at each timestep Az during the evolution is then 



AM(z) = Az / — [M(z),C,z] X 



^ M{z) 



1+e 



di (B3) 



with ( 1 + C)^(z) the mass of the (smaller) merging halo. The solid line in the left panel of Figure lBTl is the total mass accreted 
in mergers. As clear from the bottom panel, which shows the ratio between the cumulative mass grown by mergers and the total 
one derived from equation ( IBll i. mergers with 0.1 < ^ < 1 can account for most (^ 65%) of the average growth of halos. ^ 

^ We note here that the exact value of the adopted g,^i„ does not alter our overall conclusions. For example, setting ^^j,, = 0.01, increases the growth of the 
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Fig. B 1 . — Left: Growth history of halos with final mass ~ 1 — 2 X IO'^Mq, as derived from the halo mass function (long-dashed fine), compared to the 
growth history as predicted by integrating the merger rates of halos with progenitors mass ratios in the range 0. 1 < ^ < 1 {solid fine). Right: the long-dashed 
line shows the total growth in the central black hole as derived from the long-dashed line in the left panel and using equation jsj; the solid line is instead the 
black hole growth derived by integrating the central galaxy merger rates. The dot-dashed line in the lower panel con'esponds to black hole growth when no delay 
is considered. Mergers contribute by only < 7% to the overall growth of the central black hole, at the most, if the delay due to dynamical friction is taken into 
account. See text for details. 



It is then natural to ask how muc h of the central black hole mass is actually contributed by mergers. The long-dashed line in 
the upper right panel of Figure IbTI shows the total growth of the central black hole derived by assuming that at each z the black 
hole has, on average, a mass as given by equation (|5]l. In order to estimate the contribution of black hole mergers we, however, 
cannot naively use equation ( IB3l l. We need to in fact take into account that when the smaller halo enters the virial radius of the 
larger halo, it takes about a dynamical friction timescale to sink to the center allowing for the central galaxies (and their black 
holes) to actually merge. 

We therefore follow Shen (2009), and compute the central galaxy merger rates as Bgai[A^(z),^,z] = being z 

the redshift of the actual merger with the central galaxy and Ze the redshift at which the smaller halo first entered the virial radius 
of the larger halo. We use the Shen (2009) analytical fit to the function which reproduces the Jiang et al. (2008) merger 

timescales well in the range 0.1 < ^ < 1. Adopting the results by Jiang et al. (2008) is particularly meaningful for this paper 
Their subhalo merger timescales were in fact derived for a suite of high-resolution numerical simulations performed on halos 
with masses M > 5 X W^H-^Mq, in the range of interest for our paper, and with a virial mass definition equal to that used in 
this paper (see §|2]). The study by Boylan-Kolchin et al. (2008), for example, yields somewhat different merging timescales with 
respect to those found by Jiang et al. (2008; see also Shen 2009). However, their results were based on parent halos about an 
order of magnitude less massive than the ones of interest here, and with a significantly different definition of virial mass. 

The solid line in the upper right panel of Figure lsTl shows the total mass accreted onto the central black hole via mergers. It is 
clear that subhalo mergers are not the dominant source of growth for massive black holes at high redshift, as also independently 
found by several other groups (e.g., Volonteri et al. 2005). The cumulative fraction of black hole mass at z = 3 grown via mergers 
is only ^ 7%, reducing to just a few percent at 4 < z < 4.5 where most of the tightest constraints from clustering come from 
(see §[3]|. Adopting the Boylan-Kolchin et al. (2008) merging timescales would increase the contribution from mergers by just a 
factor of < 2. Note, also, that our estimate is actually an upper limit. In fact, this calculation assumes that all the incoming dark 
matter substructures actually contain a black hole as massive as what predicted from equation (|5]l that efficiently merges with 
the central black hole. Moreover, we have not considered that the delay time for black hole mergers is even longer than those of 
galaxies (see, e.g., Merritt & Milosavljevic 2005 for a review), a correction which would further drop the contribution of black 
hole mergers. 

In the lower right panel of Figure IB II we also show the predictions of the same model when no delay is considered between 
halos and central black hole mergers. In this model, satellite black holes instantly merge with the central black hole of the parent 
halo just after the merging of their host halos. Although this assumption is obviously too simplistic, as it neglects any dynamical 
friction time delay, it can be safely regarded as a secure upper limit to the contribution of mergers to the overall black hole growth. 
As shown in the right panel of Figure IbTI in this extreme model the growth of black hole mass via mergers is comparable to that 
via accretion, accounting for about 50-60% of the final black hole at z ^ 3. 

From the study undertaken in this section, we conclude that, in the physically plausible case that a significant dynamical friction 
time-delay is present between host halo and central black hole mergers, it is a good approximation to neglect black hole growth 
via mergers in the continuity equation model discussed in §|2] However, the same is not true for quasar activations. The model 
discussed in §|2]holds in reproducing both the quasar luminosity function and quasar clustering only in the hypothesis that black 
hole growth via accretion parallels that of the host dark matter halo with no time delay between the two. The latter assumption 
was also adopted by several other groups to boost and thus facilitate the match to the high-z, luminous quasar number counts 
(e.g., Wyithe & Loeb 2003; Scannapieco & Ho 2004; Lapi et al. 2006; Shen 2009). 

For completeness, however, we present in the next section the results of a fully self-consistent model that evolves the black 
hole mass function through a continuity equation with accretion and mergers, with no delay in any of its components. We will 
discuss how the outcome of such a model strengthens our general conclusions. 



parent halo via mergers to 90 — 95%, but it has a negligible impact on the mass growth of the central black hole. 
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MODELS WITH ACCRETION PLUS MERGERS 
Inserting mergers in the continuity equation of equation ( fTOl l implies a format of the type 



^»bh(Mbh,0 
dt 



1 



5$(L,z) 



fefln(10)2MBH 91ogL 



where 



Af 



,2 



merg.z 



Af 



, ^BH 



and 



Af 
p 

J me 



Mbh,z 



Af 



:(e,M)«4M(M{,'H = (i+e)MBH, 



dMl^ dMBH 



dM c/M^jH 
■J dM^dM^ 



(CI) 



(C2) 



(C3) 



are, respectively, the merger rate of incoming smaller mass black holes with mass Mg jj — MbhC / ( 1 + C) ™d Mg jj = Mbh / ( 1 + 
that merge into a black hole of final mass Mbh, and the merger rate of black holes with initial mass Mbh that merge into more 
massive black holes of mass M'^^ = Mbh(1 +0/$ ™d ■'^bh — -^bhII +0 (in both Eqs. [C2] and [C3] we set ^^in = 0.1, and 
add the factor of 1 /2 to avoid double counting). 

If we assume that no delay is present between the mergers of the black holes and their parent halos, the probability of black 
hole mergers per unit time is simply given by the halo merger rate adopted above, i.e.. 



Af 



i {i,M)nh [M(Mbh ,z),z]^Bh [M, ^(z) , z] . 



(C4) 



By simply knowing, at each timestep, the mapping between infalling halo mass and its central black hole (given by Eq. ||5l), we 
can then compute the expected average rate for any black hole merger event. By further assuming the AGN luminosity function 
to be known from observations (we here adopt the analytical derivation by SWM), we can simply integrate the right-hand side of 
equation ( ICll ) and derive the black hole function at all redshifts. 

The result is shown in Figure ICTl the left panel of which shows the integrated black hole mass density /9(Mbh > lO^M0,z), in 
the mass and redshift range of interest here, for the reference model (A = 0.25, a = 1, e = 0.15) with no mergers (long-dashed 
lines), and with mergers (solid lines), as labeled, while the right panel compares the resulting black hole mass functions for the two 
models at three different redshifts, from bottom to top, z = 6, 4, 3. The two filled circles in the left panel mark the expected black 
hole mass density at z = 3 and z = 4 expected from equation ( |20] |. adopting a duty cycle of Pq — 0.25 and Pq = 0.5, respectively, 
for a model consistent with the measured quasar clustering (see Figure [T2li. As above, we here neglect any source of scatter in the 
relation between luminosity and halo mass to maximize the predicted clustering for a given model. We find that, irrespective of 
the differences in the AGN luminosity function adopted here and in the main text, the results for the pure accretion model match 
those presented in Figures [12] When mergers are included, the estimated black hole mass density above 1O^M0 is larger by a 
factor of > 2, than the one from accretion alone. We should emphasize that we consider this factor of two to be an extreme upper 
limit on the impact of mergers, since it ignores even the effects of dynamical friction delay, which we have shown in the previous 
Appendix to drastically reduce the black hole merger growth. More importantly, however, any growth of the black hole mass 
function via mergers exacerbates the tension we have highlighted between high-redshift quasar clustering and the luminosity 
function. Mergers add mass to the black hole population without associated luminosity, so reproducing the observed luminosity 
function with the black hole population implied by the black hole-halo relation requires a still higher radiative efficiency. 
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Fig. C 1 . — Left: accreted black hole mass density above logMen/^© > 8 as a function of redshift, for the reference model characterized by A = 0.25 and 
e = 0.15 , with no mergers {long-dashed line), and with mergers {solid line); the solid, filled circles are the z = 3 and z = 4 black hole mass density obtained via 
Eq. 1201 by assuming a duty cycle of Pq = 0-25 and Pq = 0.5, respectively, as in the reference model (Figs.[8land ll2) . Right: comparison between the resulting 
black hole mass functions for the no-merger {long-dashed) and merger (solid) models, at three different redshifts, from bottom to top, z = 6, 4, 3. 
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